function p = ls_2Dint(b, f, mat,e, angt, Dt0, x, z, varargin)
% p= ls_2Dint(b, f, mat, e, angt, Dt0, x, z, Nopt)computes the normalized 
% pressure, p, for an element in a 1-D array radiating waves across 
% a plane fluid/fluid interface where p is calculated
% at a location (x, z) (in mm)in the second fluid for a 
% source of length 2b (in mm) at a frequency, f, (in MHz).
% The vector mat = [d1, c1, d2, c2] where d1 is the density in the first
% medium (in gm/cm^3), c1 is the wave speed in the first medium
% (in m/sec)and similarly d2 is the density in the second medium (in 
% gm/cm^3)and c2 is the wave speed in the second medium (in m/sec).
% The distance e (in mm) is the offset of the center of the element from
% the center of the array. The parameter angt(in degrees)
% specifies the angle of the array with respect to the x-axis 
% and Dt0 (in mm) is the distance of the center of the array from the 
% interface. The assumed harmonic time dependency is exp(-2i*pi*f*t).
% The model used is a Rayleigh-Sommerfeld type of integral for a
% piston source where ray theory has been used to propagate the cylindrical
% waves generated by the element across the interface.
% Nopt gives the number of segments to use. If Nopt is not
% specified as an input argument the function uses one segment 
% per wavelength, based on the input frequency, f, which must
% be a scalar when Nopt is not given. 

% extract material parameters
d1 =mat(1)  ;
c1 = mat(2);
d2 = mat(3)  ;
c2 = mat(4);
% compute wave numbers
k1b = 2000*pi*b*f/c1 ;  
k2b=2000*pi*b*f/c2;

% if number of segments is specified, use

if nargin == 9
    N = varargin{1};
else
% else choose number of segments so that the size of each segment
% is a wavelength 
    N = round((2000)*f*b/c1); 
    if N < 1
     N=1;
    end
end

% compute centroid locations for the segments
xc =zeros(1,N);
for jj=1:N
    xc(jj) = b*(-1 + 2*(jj-0.5)/N);
end

% calculate normalized pressure as a sum over all the segments

p=0;
for nn= 1:N
    % find the distance, xi, where the ray from the center of a segment 
    % to point(x,z)intersects the interface 
    xi = pts_2Dintf(e, xc(nn), angt, Dt0, c1,c2, x, z);
    % compute distances and angles needed in the model
    Dtn=Dt0+(e+xc(nn)).*sin(angt*pi/180);
    Dxn = x-(e+xc(nn)).*cos(angt*pi/180);
    r1 = sqrt(xi.^2.+ Dtn.^2)./b;
    r2 = sqrt((Dxn -xi).^2 +z.^2)./b;
    ang1 = asin(xi./(b*r1));
    ang2 =asin((Dxn-xi)./(b*r2));
    ang = angt*pi/180 -ang1;
    ang = ang + eps.*( ang == 0);
    % form up the segment directivity
    dir =sin(k1b.*sin(ang)/N)./(k1b.*sin(ang)/N);
    % compute plane wave transmission coefficient(based on pressure ratio)
    Tp = 2*d2*c2.*cos(ang1)./(d1.*c1.*cos(ang2) +d2.*c2.*cos(ang1));
    % compute phase term and denominator
    ph =exp(1i*k1b.*r1 + 1i*k2b.*r2);
    den =r1+(c2/c1).*r2.*((cos(ang1)).^2)./(cos(ang2)).^2;
    % put terms together for pressure due to each segment
    p= p + Tp.*dir.*ph./sqrt(den);
   
end
    p   = p.*(sqrt(2*k1b./(1i*pi)))/N;  % include external factor

